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SUMMARY 


Nonlinear beam kinematics are developed and applied to the dynamic analysis of 
a pretwisted, rotating beam element. The common practice of assuming moderate rota- 
tions caused by structural deformation in geometric nonlinear analyses of rotating 
beams has been abandoned in the present analysis. The kinematic relations that 
describe the orientation of the cross section during deformation are simplified by 
systematically ignoring the extensional strain compared to unity in those relations. 
Open cross section effects such as warping rigidity and dynamics are ignored, but 
other influences of warp are retained. The beam cross section is not allowed to 
deform in its own plane. Various means of implementation are discussed, including a 
finite element formulation. Numerical results obtained for nonlinear static problems 
show remarkable agreement with experiment. 


1. INTRODUCTION 


It is now widely recognized that aeroelastic analysis of helicopter rotor blades, 
particularly of hingeless and bearingless rotor blades, requires the incorporation of 
kinematical nonlinearity (ref. 1). The main reason for this requirement is that the 
stability and response of such systems depend strongly, in some cases, on the cou- 
pling between bending and torsion motion. This important coupling cannot be obtained 
accurately without consideration of the kinematical nonlinearities. 

A brief history of the developments of nonlinear equations of motion for rotat- 
ing beams prior to 1974 is given in reference 1. These developments are generally 
concerned with slender beams, with the effects of shear deformation ignored. Since 
1974, the major contributions to this subject have been those of Kaza and Kvaternik 
(ref. 2), and Rosen and Friedmann (ref. 3). The equations of motion developed in 
references 1-3 are very similar. In references 1 and 3 an ordering scheme is used 
that limits the kinematical development to moderate rotations. In reference 2 the 
nonlinearities are limited to the second degree in the displacement variables. These 
two different methods of specifying "moderate rotations" in references 1-3 are vir- 
tually equivalent for this problem. 

In a general-purpose analysis, a single set of equations is desirable - one that 
is valid for all values of the equation parameters, within some range. When moderate 
rotations are assumed, situations can easily arise in which the solution violates the 
assumption of moderate rotations. One example is the case of a thin beam for which 


the ratio of bending stiffnesses is small compared to unity. In order to avoid this 
problem, certain ad hoc modifications to the equations of reference 3 were intro- 
duced and were necessary in order to produce the excellent correlation obtained in 
reference 4 with experimental data for large displacements of an end-loaded canti- 
lever. For example, magnitudes of the beam bending and torsion stiffnesses had to 
be specified before the equations could be put into final form for solution. 

Ideally, the magnitude of parameters in the equations for general-purpose analyses 
should not influence the equations themselves. Such an ideal is evidently not pres- 
ent in the ordering schemes of references 1 and 3 or in any arbitrary a priori 
restriction to second-degree nonlinearity, as in reference 2. 

There are other shortcomings in the equations in references 1-3. For example, 
the effects of pretwist are not treated rigorously. An improved treatment of pre- 
twist effects is presented in reference 5 for a simplified problem involving only 
torsion and axial displacement. Additional work is required to incorporate those 
analysis techniques into a general, nonlinear, bending-torsion-extension analysis. 

In reference 6 exact nonlinear kinematical relationships are developed and additional 
insight is presented concerning relationships among the equations of references 1-3. 
Finally, in reference 7 it is shown that inconsistencies are virtually unavoidable 
in ordering schemes based on displacements and rotations when the magnitude of the 
torsion rigidity is small compared to bending stiffnesses. 

Because of the problems with kinematical limitations in the above approaches, it 
seems appropriate to model the kinematics of a slender beam without resorting to an 
ordering scheme on rotations or to arbitrary restrictions on degree of nonlinearity 
allowed in expressions involving displacement. This method would avoid some of the 
limitations of previous analyses and circumvent nonrigorous modifications such as 
were found necessary in reference 4. In the present work, the development in refer- 
ence 6 serves as a foundation along with some important observations from references 8 
and 9. Rather than develop the partial differential equations of motion for this 
problem, the objective is to develop a statement of the principle of virtual work for 
dynamic analysis of a rotating beam element with Euler-Bernoulli kinematics. This 
statement will serve as the basis of a Ritz-type modal model for an entire rotor blade 
or of a single, finite element of a rotor blade without having to consider the partial 
differential equations and the natural boundary conditions, which are available from 
the statement, if needed. 

Although the results presented here are for the case in which shear deformation 
is neglected, the necessary relations for including shear deformation are included as 
part of the development. Similarly, with the appropriate constitutive law, effects 
such as orthotropy or anisotropy could be incorporated. The initial curvature of the 
elastic axis and effects associated with open cross sections could also be incorpor- 
ated. The detailed development of these topics is reserved for future extensions of 
the present analysis. The present development proceeds as follows: The beam kine- 

matics are developed in section 2 in such a way that the necessity for an ordering 
scheme is obviated. This section includes the development of strain-displacement 
relations and direction cosines of the nominal cross-sectional plane. In section 3 
generalized forces caused by internal loads are developed from a consideration of 
the strain energy. Expressions for generalized forces caused by inertial and gravi- 
tational loads are developed in section 4 based on the work done by these loads 
through a virtual displacement. Generalized forces caused by a general set ol 
applied distributed loads are developed in section 5 similar to the development in 
section 4. A finite element implementation is discussed in section 6. Finally, 
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numerical results obtained from a finite element calculation of nonlinear static 
equilibrium are presented in section 7. 


2 . KINEMATICAL DEVELOPMENT 


In this section, the kinematics for the beam element are developed starting with 
the rigid-body motion of the cross-section plane as described in reference 6. Cross- 
section warp is then superimposed on the rigid-body motion to obtain the final dis- 
placement field. Next, the strain-displacement relations are developed from the 
displacement field. The extensional strain is then assumed to be small compared to 
unity. This assumption is used to simplify the orientation description and moment 
strains considerably, without sacrificing accuracy. 


2.1 Development of Displacement Field 

First, consider a straight-beam element with the associated coordinate systems 
shown in figure 1. It is assumed that the motion of the frame F is known in an 
inertial frame I. A set of dextral axes x-^, i = 1,2,3 is assumed to originate at 
F*, the origin of F. The x 3 -axis lies along the elastic axis of the beam element. 
Each unit vector b?, i = 1,2,3 is parallel to the corresponding axis xi. The beam 
is assumed to be pretwisted so that the local-beam cross section is rotated by an 
angle 0 about the x 3 -axis at any point on the x 3 -axis. At x 3 = 0, 0 is defined 
to be zero so that the major axis is parallel to the x 2 -axis and the minor axis is 
parallel to the x^axis. Consider an arbitrary material point in the beam prior 
to deformation, denoted by M c . The position vector of M 0 with respect to F x is 

MnF* F F F 

R 0 = x 3 b 3 + (5 1 cos 0-5 2 sin 0)b 1 4- (E >1 sin 0 + K 2 cos 0)b 2 (1) 

Denote the same material point in the deformed beam by M. The position vector of 
M with respect to F* is 

MF* F F P P 

R = x b + u.b. + 5 b + ¥(x JAU ,£jb a (2) 

- 3-3 1-1 a-a 3 12 — 3 ' 


where repeated Latin indices imply summation from 1 to 3 and repeated Greek indices 
imply summation from 1 to 2 . 


The axes E >1 and £ 2 are along principal axes for the local cross section at a 
point P* on the elastic axis and remain so during deformation. The unit vectors 
are parallel to the cross-section principal axes at P*, which islocated at the 




origin of a dextral system P. The unit vector b 3 is defined as b^ x b 2 and is 
thus normal to the cross section. The displacements u^, i = 1,2,3, are along 


bj. 


respectively. A warp displacement field has been added vectorially to the rigid-body 
component of the position following Wempner (ref. 9). The warp amplitude is 4 , and 
X is the cross-section warp function. 


The basis vectors for the P axes are denoted in equation (1) by 


P F 

b. = C. .b. 
-i 1J-J 


(3) 
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where Cj_j is a 3 x 3 matrix whose elements may be specified by any of several sets 
of parameters (ref. 10). All three-parameter descriptions of the rotation have 
inherent singularities. Classical Euler orientation angles (called body-two orien- 
tation angles in reference 10) have singularities at values of certain angles equal 
to zero. The Tait-Bryan orientation angles (called body-three orientation angles in 
reference 10) have all singularities at 90°. This fact makes them more amenable to 
descriptions of rotation caused by structural deformation than Euler angles, since 
the case of no rotation would correspond to no deformation and Euler angles are 
singular at that condition (ref. 6). To avoid singularities altogether, at least 
four parameters are required to describe the rotation. One of the best known 
descriptions is the set of Euler parameters. Euler parameters enable rotation to 
be described with four scalar quantities upon which the direction cosine matrix ele- 
ments depend quadrat ically . It is possible to eliminate one of the Euler parameters 
algebraically and derive the Rodrigues parameters. There is a singularity at 180° 
of rotation in any direction, but this is rarely a problem in deformable structures. 
The direction cosines are simple ratios of quadratic polynomials in the Rodrigues 
parameters. Furthermore, the inverse operation (i.e., given the direction cosine 
matrix, find the parameters) is trivial compared to the same operation with orienta- 
tion angles (ref. 10). 


In this paper, both Tait-Bryan orientation angles and Rodrigues parameters are 
considered. The primary development is executed with the orientation angles because 
of the simplicity of the result and the familiarity of the method. Rotations free 
of singularity up to 90° are completely acceptable in helicopter blade applications. 
For applications in which the 90° restriction is unacceptable, the kinematical 
development for Rodrigues parameters is given in the appendix. 


The direction cosines are here expressed in terms of Tait-Bryan orientation 
angles 0£, i = 1,2,3 defined as follows. Let basis vectors b? be introduced 
beginning with b? aligned with b^. Then perform sequential rotations 0i about 
b^, i = 1,2,3 until b? is aligned with the principal axes of the deformed beam 

i = 1,2,3, as shown in figure 2. Other sequences of rotations are possible, but 
here is used. Direction cosines of the local, def ormed-beam, cross- 

section principal axes P with respect to F when expressed in matrix form are 


C 


C 2 C 3 S 3 C 1 + S 1 S 2 C 3 S 3 S 1 " C 1 S 2 C 3 

-c 2 s 3 c 3Cl - Sl s 2 s 3 c 3Si + C;l s 2 S 3 


- S 1 C 2 


C 1 C 2 


(4) 


where c^ = cos 0p and s± = sin 0^. The displacement field is now completely spe- 
cified although, as pointed out in reference 6, two of the three angles and 0 2 

can be eliminated if shear deformation is neglected. 


2.2 Development of Strain-Displacement Relations 

Now that the displacement field is determined in terms of Uj[ and 0-^, it is a 
straightforward matter to calculate strain-displacement relations. There are several 
strain measures that could be used in this type of analysis. Almansi strain was used 
by Hodges and Dowell (ref. 1) in their preliminary development, and Green strain has 
been used by almost everyone else. For infinitesimal strains, however, differences 
among the common definitions of strain are small and are ignored in this 
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development. The intent here is to develop a set of strain-displacement relations 
that will be linear in elongations and shears, but unrestricted in rotations up to 
changes in orientation where singularities are encountered. Thus, it is immaterial 
whether it is Green strain or Almansi strain that serves as the starting point since 
the relations will be simplified for small strains anyway. 

M F* MF* 

The independent variables of the vectors R 0 and R (namely, , £ 2 , and x 3 ) 
constitute a nonorthogonal curvilinear coordinate system identical to the one used in 
reference 5- The steps followed in reference 5 to derive a set of strain displace- 
ment relations from these vectors are as follows: 

1. Obtain the covariant base vectors for the undeformed state from equation (1) . 

2. Obtain the covariant metric tensor for the undeformed state from step 1. 

3. Obtain the contravariant metric tensor for the undeformed state from step 2. 

4. Obtain the relationship between a local Cartesian coordinate system and 
step 1 from step 3. 

5. Obtain the covariant base vectors for the deformed state from equation (2). 

6. Obtain the covariant metric tensor for the deformed state from step 5. 

7. Obtain the Green strain tensor from half the difference between step 6 and 
step 2. 

8. Transform the Green strain tensor from step 7 to the local Cartesian coordi- 
nate system using step 4. 

9. ignore elongations and shears with respect to unity in all strain components 
(i.e., discard all squares and products of elongations and shears, leaving each 
strain component to be a linear combination of elongations and shears) . 


The details of the algebra, although lengthy, are straightforward and are 
omitted in this report. The engineering strain components resulting from these 
operations are 


C 31 = e 3 


! + (X x - ? 2 ) (< 3 - 6 ') 

E 32 ~ E 32 (^2 ^ 1^ K 3 - ® ) 




> 


(5) 


M ■ Ss + Mi - Mi + M + MM - e "> 12 

+ (5 z x 1 - SjXjJCk, - e ’ ) a ' + x(k 3 - e ■) ’ 

The strains at the reference axis, referred to as force strains, are given by 

c 3 i = C-j_j(6 3 j + u j) “ ^ 3 i (6) 


where 6ij is the Kronecker symbol. The curvature-like quantities, referred to as 
moment strains, are given by 
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(7) 


"i ‘ “U 9 i 


where, in matrix form. 


(JO 


-c 2 s 3 C 3 0 
. S 2 01 


( 8 ) 


Derivatives of the warp function are denoted by A a = 3A/3£ a . A restricted warp 
amplitude I = K 3 - 0 f is assumed here instead of the more general formulation in 
terms of I explicitly as in reference 5. This assumption results in the neglect 
of transverse shear in the outer fibers of the beam. It is interesting to note that 
the force and moment strains are very similar to those of Reissner (ref. 8 ) except 
that the ones above are expressed in terms of Tait-Bryan orientation angles instead 
of Rodrigues parameters as in reference 8 . The moment strains have the dimensions 
of curvature, but differ from curvature by a factor of s ! , where s is the length 
coordinate along the deformed-beam elastic axis (ref. 6 ) and ( ) 1 denotes the 
derivative with respect to X3. The moment strains also closely resemble those of 
Wempner (ref. 9 ) developed for small deformation of arbitrarily curved beams. 

It is possible to eliminate 0 X and G 2 from the analysis if shear deformation 
is ignored, resulting in an Euler-Bernoulli beam model. For this case the vector 
tangent to the beam elastic axis, 3 R MFx / 3 s|r _ 0 , must remain normal to the local, beam 

cross section during deformation. Thus, from equation (2) 


MF* 1 F P 

SR /Ss| = [6 3 i (9x 3 /Ss) + Su ± /Ss]b i = b 3 

By virtue of equations (3) and (9) 

3x 3 

C 3i = 6 3i ~Ys~ + 'YT 


(9) 


( 10 ) 


Since Cpj is orthonormal, C 3 ^C 3 j_ = 1. If all derivatives are expressed with respect 
to x 3 instead of s, an expression for s f is obtained: 

s' = [up + u ’ 2 + (1 + u ') 2 ] 1 / 2 ( 11 ) 

From equations (10) and (4) 

S 2 = u i y/s? ( 12 a) 

~*s 1 c 2 = u^/s 1 ( 12 b) 

c 1 c 2 = (1 + u’)/s ? ( 12 c) 

The angles 0 T and 0 ? can be eliminated from C-^j with the following relations 
obtained from equations ( 12 ) 
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(13a) 


s, = 


_u 2 


(s' 2 - u’ 2 ) 1 / 2 


1 + 


I 2 \ l/2 


_ (s' 2 - u[ z - u’ 2 ) 
(s' 2 - up) 1 / 2 (s' 2 - uj 2 ) 1 / 2 


(13b) 


s n = 


(13c) 


C 2 


(s' 2 - u* 2 ) 1 / 2 


(13d) 


Equations (12), when substituted into the expressions for force strains, yields 

z =0 
3a 

c_ = s 1 - 1 J 
3 3 


(14) 


Considerably more algebra is required to eliminate 0 1 and 0 2 from the moment strains 
if the above method is used. It is relatively simple, however, to write them directly 
from reference 6 in terms of derivatives with respect to s where ( ) + = 3/9s( ). 


(1 _ „+ 2 W 2 
- „+ 2 W 2 " (1 _ 


"l S 3 


(1 - u+ 2 ) 


[-H- -t-frr i 
1 - U«J 


(15a) 


u c 
1 3 


(i 


(l-ut 2 ) l/2 (1- 


-u +2 ) 1 / 2 s r++ U+U+U++1 

— ■ - 1 ■■ — 3 r 

u + 2 _ u +2)l/2 2 i - U+ 2 


_JL = n+ 

s’ 3 (]. - 


■f r + + -H-n 

U 1 4+ , U 1 U 2 U 1 

u+ 2 - uf) 1 / 2 p 1 - ut 2 J 


(15b) 


(15c) 


Equations (15) may be expressed in terms of ( ) f quantities instead of ( ) + quanti- 
ties by equation (11) and lengthy algebra. This is unnecessary, however, in light of 
the simplifications in the next section. 


2.3 Small-Strain Simplification of Kinematics 

Because the moment strains are complicated, it is useful to simplify them 
through the derivation of a small-strain approximation. This is accomplished by 
neglecting the longitudinal strain of the elastic axis with respect to unity in the 
direction cosines as well as the moment strains. The expression for the longitudinal 
force strain s’ - 1 is already linear in the elongation of the elastic axis; the 
moment strains that are to be obtained are then independent of elongation of the 
elastic axis (i.e., independent of s T ) . It is interesting to note that Reissner’s 
moment strains differ from curvatures by a factor of s 1 and are independent of s f . 
When 0 X and 0 2 are eliminated, additional s 1 terms are introduced into the 
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equations through equations (13). Simply setting s T = 1 in all places that it 
occurs in the moment strains and direction cosines removes the dependence of these 
quantities on e 33 . 

Another way to view this approximation is to expand the strain in a Taylor 
series with elongations caused by stretching of the elastic axis as the small param- 
eter. This separates elongation into components on and off the elastic axis. When 
only the terms linear in elongations are retained, the result is equivalent to having 
set s T = 1 in all quantities except £ 33 * Thus strain is considered small compared 
to unity as in the example given in reference 6 , p. 32. It is claimed in reference 11 
that this sort of approximation invalidates the strain for cases other than inexten- 
sional. In reality, a simplification of this sort cannot significantly affect the 
accuracy of the mathematical model as long as the strains are small relative to unity, 
which they must be for applications of Hooke's law. The resulting simplifications in 
the derivation and in the final equations are substantial. 


If s T is set equal to 1 in equations (13) and (15), the result is 



(16a) 


(16b) 


(16c) 


(17a) 


(17b) 


u. 


(17c) 


c 


2 



1/2 


(17d) 


Note that the sign of the square root quantities becomes ambiguous when 0 1 or 0 ^ 
exceeds 90°. Thus, it is imperative in this formulation to restrict orientation 
angles to less than 90 . Also note that k± and Cj_j are now independent of u'; 
s' must not be set equal to unity in the force strain c 33 = s f - 1 . 

Geometric boundary conditions are determined from specified u^ and Cpj at 
x 3 = 0 and x 3 = £. Natural boundary conditions may be obtained from the complete 
statement of the principle of virtual work as developed in the next three sections. 
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3. DEVELOPMENT OF GENERALIZED FORCES CAUSED BY INTERNAL LOADS FROM STRAIN ENERGY 


In this section, the internal loads for the beam element are developed from the 
strain-displacement relations from section 2 and a conventional constitutive law. 

The beam is assumed to be of such a configuration that open-cross-section effects, 
such as warping rigidity , are negligible, which is justified for rotor-blade cross 
sections. However, the effects of warping are retained in other parts of the inter- 
nal loads development. This assumption is helpful in a finite element context since 
enforcing kinematical boundary conditions on the warp displacement field at finite 
element nodes is not possible in the general case of beams being joined at arbitrary 
positions and orientations with respect to one another. 

The virtual work on the internal loads is obtained from the variation of the 
strain energy U as in (ref. 9) 


5U 



(Ec 6e + Ge 6e )d£ 
v 33 33 3a 3 a 1 


d£ dx 
2 3 


(18) 


where 



6£ 33 = 

<5s ’ 

+ 5 2 « k i - 5l « s + ml + e|)(K 3 

6£ 31 = 

(A i 

5 2 )<s< 3 

6£ 32 = 

( A 2 

+ 5 x ) 5 k 3 


- o’) + U x - 5 i )e f ]<$K + xskO 

2 1 12 3 3 


> 




(19) 


The expressions are greatly simplified if the variations of force and moment strains 
are written as 


x , 3s' , 

6s = 3^ 6u i 


(20) 


3k . 5k. 

- 3^ K + 4 i3 #e J + 3^ K + e ia, K a 56 , 

a a 


where e ijk is the Levi-Cevita permutation symbol and 


5s ' 
3u. 


6 3i + u! 




(21) 


8K a 

3u" 


3k 


a2 


* 3 3 


5k 

3 

Su’i' 


-C 2 C 

3132 


C (1 - C 2 ) 

33 V 3 ± 


5k . 


-C . 

:u 

C 

33 


y 


u" r C,,C,, C C 1 u"C C C 2 

1_ I 22 32 + 12 3 1 _ 2 11 31 32 

C 2 I 1 - C 2 C 3 3 I C 3 (1 - C 2 ' 

33 L 31 J 33^ 31' 


J 


( 22 ) 
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u u C C C C 

1_ 12 3 2 22 3 1 

P 2 1 n 2 C,, 

C 3 3 L 1 " C 3 1 33 J 


u"C c c;„ 

2 21 31 32 

: 2 d _ c 2 ) 

33'“ U 31 ; 


_U 1^3 1^32 ^^3 3 + ^31 C 31^ U 2^ ^32' 

3 2 2 3 

C (1 - c ) C 

33 31 33 


3k . C 
1 11 


( u l^ 31 U 2^ 3 2^ 


Next, the strains are substituted, equations (5) with e 3 j_ given by equation (14), 
k i by equations (16), and the virtual strains given in equations (19) into the 
strain energy by equations (18). The resulting expression can be arranged by 
terms that multiply 6s 1 , 6k^, and 6k’. These coefficients are denoted by F 3 , Mj_, 
and B, respectively, which are given by 


F 3 = E 0 (s’ - 1) + E 2 k x - E x k 2 + -J- ( k 3 - 6’) 2 + D o 0'(k 3 - 6’) 


6 0 I d k 0 + e D „ 
ag g 3 af$3 




1') + ~ (K 3 - e’) 2 + Dg6 ' (k 3 


3D,e ' 


- e')J 


J + I 3 < 8 ' - 1} - e aP,Vp + T (K 3 " 01)2 + — 2 ~ (K 3 - 0,) + D 4 0 ' 2 


x (k 3 - O') + [D 0 (s- - 1) - e a33 D a K 3 ]0’ 


B = S (k - 0') + S 0' 

13 2 


These quantities are the stress resultants given in terms of the displacements. The 
distributed axial force is F 3 ; local distributed moments along b? are denoted by 
Mj[. The distributed bimoment is B, which may be neglected for the class of beams 
considered herein. The section properties used in equation (23) are defined by the 
following integrals over the cross section. 


JI e «> 


D q “ II E ^2 


J a E ^a ^i ^2 ’ - JT E ^a e $y3^B^Y °^2 

A •'•'A 

: 1 2 JI E{ 3 dt . «» ; °3 E JI E5 yS 'WaS «, d «: 
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I 2 = JJ E5 2 d? x d£. 


= // d5 i d ^ 


B a = IT d5 x d? 2 ; J 5 JT G[(^ + X 2 ) 2 + U 2 - d5 s 

"A "A 

B , 5 JT E <W 2 d «i d ?2 i JJ E Ma d «l «2 s 0 i b ■ *1 + *2 

A A 

s 2 - JJ a ex <Mi - «i X 2> d «i d 5 2 

s x = JJ EX£ a £ a ^1 ^ 2 

JJ E\E, a dq d^ 2 = 0 ; JJ EX d£ 2 E 0 

The strain energy is then 

* r h + Mi ^t) k + f > ^ ^ + 


(24) 

(cont) 


-™i„t - 5D - 


8 k . -| 

+ M. T-n- <5u" + M 3663 dx, 

1 dU a 6 0 I 4 

a J 


(25) 


A suitable discretization of and 0 3 is sufficient to obtain generalized forces 

in a Ritz-type formulation. 


4. GENERALIZED FORCES CAUSED BY INERTIAL AND GRAVITATIONAL LOADS 


In this section the generalized forces caused by inertial and gravitational 
loads (i.e . , body forces) are developed for the beam element. The effects of cross- 
section warp on the body forces are negligible and are not considered. It is 
assumed that the gravity field g, the velocity of point F* in an inertial frame 
y^* 1 , the acceleration of point ~F* in an inertial frame a F *I, the angular velocity 
and angular acceleration of frame F in an inertial frame and respectively, 


FI 


in frame 

F as 



S F i-i ’ 

F*I 

a 

F*I, F 

= a b . ; 

Fi -1 

F*I F*I, F 

v = v_ . b . 

Fx -1 

FI, F 

: -Fi-i 5 

FI 

a 

FI F 
= a„ # b . 
Fi-i 



(26) 
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We must find the acceleration and virtual displacement of the point M where the 
position of M with respect to F* is governed by 


MF* 

R = ^3 + u + § 


where 


2s 




u.b F y 

i-i ' 


K = 


5 b P 
a-a 


( 27 ) 


(28) 


The velocity of M in an inertial frame is then 


MI F*I , FI . , N , F. , PI „ , P; 

v = v + m x (x +u)+ u + w X C + K 

TT • P * 

where r ( ) is a time derivative in F. Cross-section distortion £ is to be 
neglected. The angular velocity can be written as 


(29) 


PI PF , FI 

W =0) +03 


(30) 


Differentiation yields the angular acceleration 


PI FI FI PF , PF 

a = a +o) xo) +a 


(31) 


Now the acceleration can be written in terms of previously developed expressions as 

MI F*I FI FI FT FT F F.. PI 

a = a + a x (x g + u) + m x x (x 3 -h ia) ] -h 2 oj x u + u + a 

+ a PI x £ + i/ 1 x (w PI x p (32) 

Substitution of equations (30) and (31) into (32) yields 

MI F*I FI , , FI , FI , , , ^ o F! v F. F** 

a = a + a x (x 3 + u + + a) x [ qj x (x 3 + u + £) ] + 2w x u + u 

, 0 FI / PF . PF PF ^ , PF , FI , PF r 

+ 2 o) x (aj x p + a) x ( w x ^) + a) x (ca x Q + a x £ (33) 

The virtual displacement is obtained by the replacement of ^(*) with ^6 ( ), a vari- 
ation in the frame F, in equation (29) and the substitution of _6R and Sip for 
v and 03, respectively, in equations (29) and (30), as in reference 6, which yields 


MI F*I 

6R = SR + 


Sip 


FI 


(53 + H + 


?) 


+ F 6u 


PF 

+ Sip x Z 


(34) 


where 6R F * 1 = 5Rp^b| 
motion of the frame in 
the cross-section axes 


and Sj/ 1 = 6i(3p^b|. These quantities reflect any nonprescr ibed 
inertial space. The angular velocity and virtual rotation of 
P in F (ref. 6) can be written as 
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( 35 ) 


PF 




and 


6^ 


PF 


/3k. \ 

= 6u a + 6 3i 66 3j^i 


(36) 


The angular acceleration of P in F, a PF , is simply the time derivative in F of 
, .PF 


PF F PF /^ K i .. •• 3 k. \ p 

5 ■ * - u “ + + KH “ V-i 


The generalized forces are now found from the virtual work 

>£ 

'o J'A 

»£ 

- I 


~ 6W body " / JT P( ^ MI ‘ 8 ) ' ^ MI d5 i d? 2 dx s 


J JJ P((a MI - g) • (SR ^* 1 + F 6u) + [(x 3 + u + O 


(37) 


MI FT MT PF 

x (a - g)] • 6jj_ + [§ x (a - g) ] * 6\p Jd^ d£ 2 dx c 


(38) 


The generalized forces associated with frame motion are useful in multibody/finite 
element applications in which it is necessary to couple bodies together that are 
moving relative to each other. These generalized forces are the coefficients of 

and and can be obtained in a straightforward manner; the calculation of 

these forces is left to the reader. For the generalized forces associated with 
6uj_, 5u^, and S0 3 , the necessary ingredients are suitable discretizations for u^ 
and 0 3 and the virtual work, which is given by 


-6W, 


body 


* 5 

V r / _ F*T FI FI FT FT 

b i * J $u.^m{a - g + a x (x + u) + m x [aj x (x + u) ] + 2m 


X F U + F ii} + a FI X + U FI x (u, FI x m ? ) + 2 uj FI • m r (/ F - 2u> FI • a) PF 

. PF , PF . PF . , PF FI PF FI \ 

+ ot x nr H- co x (m x m ) + m • m m — m • m m \dx 

~ ~ / 3 


■5 

-I ' f o 6u “ + 6 3 i 60 3)^5 x k 


F*I , FI .FI 

- g + a x (x 3 + u) + a) 


FI 


x [u) J " L x (x 3 + u) ] + 2 oj F ^ x F u + F ii} + - a FF • i + (uj r ' L • i) x to 

, r>, FI .. PF . , PF . N FI , . PF . . PF . PF PF \ 

+ 2(0) • 1) x m + (m • 1) x m + (m • 1) x m + 130* -a • i\dx3 
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FI 


FI 


where the section integrals are defined as 


m 


i 3 5 jj p| • 5 d ?1 d { 2 


and 


5 JJ A 0 d L • b s K 

= JJ p§ dC, d? 2 ; i = JJ p5§ 


F. . F F.. .. F 

u = u . b . ; u = u.b. 

- 1-1 - i-i 






J 


(40) 


(41) 


5. DISTRIBUTED APPLIED LOADS 


In this section, generalized forces caused by distributed applied loads are 
simply stated for completeness. The applied force F and moment M act on the 
elastic axis so that the virtual work is 


F*I . F ( 


.FI , o , PF >. 


-6W = - f {F • (SR 1 ”' 1 ’ + A (Su) + M • (6i|/ " + A ) 


applied loads 


FI 

+ [(x 3 + u) X F] • }dx 3 


(42) 


or 


-6W 


applied loads 


= -/ m • ( 9 ^ 5u i + 6 3i 6 e 3Vi dx 3 - / e • F An 

q \ a / 0 


dx„ 


F ^ I 

- 6R 


•jL 

m n 


dx - 5^ 


FI 


r 


[M + (x 3 + u) x F] dx 3 (43) 


As with the inertial loads, the generalized forces can be obtained from the virtual 
work and any suitable discretization of u^ and 0 3 . 

6. VARIOUS SCHEMES OF IMPLEMENTATION 


The total virtual work may be written as 


<SU 


6W - 6W . . , , , = 0 

body applied loads 


(44) 


where the various components of the total virtual work may be obtained from equa- 
tions (25), (30), and (43). The kinematical variables u-^, i = 1,2,3, and 0 3 may 
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be discretized in several ways. One way is to use a set of admissible functions in 
a Ritz-type analysis based on equation (44) where all the frame variables are pre- 
scribed and 6RF*1 and are zero. A variation on that type of development would 

be to allow ~6 rF*I and to be nonzero, thus collecting frame (i.e., rigid body) 

forces and moments from equations (39) and (43) directly. In this case there would 
be prescribed and nonprescribed components of each of the kinematical quantities in 
equation (26). There are clearly other ways of using the frame motion to advantage. 

If a finite-element discretization is used, it should be based on variable-order 
shape functions (refs. 12 and 13) to avoid poorly approximating the geometric stiff- 
ening term which, in this analysis, is calculated from the longitudinal strain 
s’ - 1. If u 3 is crudely approximated, this term will be inaccurate. It should 
be noted that in a redundant structure undergoing finite deformation, the geometric 
stiffening effect must be calculated from the strain-displacement relations. 

In a finite-element implementation the frame motion and forces may be used to 
advantage when coupling elements together that are defined in different moving coor- 
dinate systems, such as at the interface between rotating and nonrotating components 
of a helicopter or at a hinge. 


In figure 3, the beam element is shown with the frame F and two nodes R and T 
at the root and tip of the beam, respectively. The displacement and orientation of 
the beam cross section as a function of x 3 is represented by the displacement and 
rotation of the nodes R and T and by a variable number of generalized coordinates 
that are kinematically uncoupled from nodal translations and rotations. The beam 
displacements at the root and tip, then, must be determined from the nodal 
displacements. 


The bending displacements u 1 and u 2 are to be expanded in C x -type functions 
and u 3 and 0 3 in C°-type functions gj_, namely 


N a +i 


u = V q . (t)T . (x) 
a ai x 


i-i 


N 3 +i 


U 3 = £ ^3i (t)e i (x) 


> 


1=1 


N 4 +i 


e 3 = E %i (t >M x > 

i=i 


(45) 


where x = x 3 /i, 0 < x < 1, and Nj_ and are the orders of the shape function 

polynomials, N a > 3 and N 3 ,N 4 > 1. The generalized coordinates qjj_ with 
j = 1,2,3 have the dimensions of length and q 4 ^ is dimensionless. The boundary 
conditions on the functions are 


Vo.t) = q al (0 ; 


u 3 (0>t) = q 3i (t) ; 

u 3 U,t) = q 32 (t) 

e 3 (o,t) = q 41 (t) ; 

6 (£,t) = q (t) 
3 4 2 


( 46 ) 


15 


and on the first derivatives 


u a<0 ’ t> = V (t)/£ ’ u a (£,t) = q a4 (t ' ) ' /£ 


(47) 


The standard linear displacement functions satisfy these conditions for g: 

B, = 1 - x 


= x 


and the standard cubic displacement functions satisfy these conditions for T: 


A 

= 1 - 

3x 2 + 

2x 


= X - 

2x 2 4- 

2x 

*3 

= 3x 2 

- 2x 3 



= -x 2 

+ x 3 



(49) 


If N a exceeds 3 or N 3 ,N 4 exceed 1, the following higher-order shape functions 
allow extra generalized coordinates to be introduced (ref, 13) but still to fulfill 
the above end conditions 


3 


n+3 


= x ( 1 - x)G (5,3,x) 


'i'n-f 5 = x 2 (l - x) 2 G n (9,5,x) 


n = 0,1, . . . 


(50) 


where G n is a Jacobi polynomial (ref. 14). At the ends, the displacement and rota- 


tion of the 

nodes are related to generalized 

coordinates from 

simple kinematics. 

Let 


R R R 

a ■ u Ri-i 

5 

T 

u 

= u 

T T 
Ti-i 


(51) 


b R . b F 
-x -x 

> 

b T 

-X 

= C 

TR, R 
. . b . 


(52) 

where C TR 

depends on the pretwist 0(£) of 

the 

beam element 

at the tip of the 


element . 


c e 

s e 

o' 





c TR = 

_s e 

c e 

0 



(53) 



_ 0 

0 

1_ 





where sq = sin 0(£) and cq = cos 0(£). At the root of the element, the displace- 
ments of the node and beam are identical 


R , R R 

u b. = q. b. = q. b. 
Rx-i ^n-i M n-i 


(54) 


or 


R 

. = q . 
Rx M xi 


( 55 ) 
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Similarly, the rotations must be the same The direction cosines of R T with 
respect to its undeformed position R, R are 


c R,R (t) = C(0,t) 


( 56 ) 


At the tip, the displacements of the node and beam are identical 

T T F F F 

u b . = q, b + q b + q b 
Tl-l 1 3-1 n 23-2 32-3 

or 


(57) 


iji 

u m 

Ti 

T 

u m 

T2 

T 


r q ^ 

H 13 


V = c TR < 


q 23 ' 


q 

032 >J 


(58) 


and the rotations of the tip node are identical to those of the beam tip 

CU.t) 


C T ’ R = C T V R = C(U) 


(59) 


Equations (56) and (59) each represent nine equations, but only three are independent, 
Three preferred elements to equate are C 21 , C 3l , and C 32 since 


c 21 = -(1 - u’ l 2 ) 
C 
C 


2\ l/2 




sin 0 r 


'31 = U i 


'32 = U 2 


J 


(60) 


and are thus easily expressed in terms of q's at the root and tip. The use of 
equations (46), (47), and (60) yields for the root, from equation (56) 


l a2 _ R t R 

£ " L 3a 




q 41 = "Sin 


-l 


.R’R 

21 


(i - c?;**) 


2 \ 1 / 2 


(61) 


J 


and for the tip, from equation (59) 
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_ c t;r 2 ) 2 ' 2 


The forces and moments at the root and tip nodes must be determined from the 
generalized forces at the beam root and tip, respectively. First note that the vir 
tual displacement and virtual rotation of the root node must be identical to those 
of the beam root 


6u (t) = 5u(0,t) 
6ip R (t) = 6iK0> t) 


where 6\p = 6ip_ (see eq . 36). Similarly at the tip 


6u ( t) = 6u(Z, t) 
( t ) = 6ip(£,t) 


The relations for the virtual displacements are simple first variations of equa- 
tions (35) and (58). The rotations are more involved. The use of equation (36) 
evaluated at the root yields 


(0 ’ 


~ (0,t) — P- + <5 6 q 


si l ij 


C. .(0,t) 


where . In matrix form 


R(0,t) J Sip* 


where 


R(x 3 ,t) = 


-£C,., 


1 - Cfi 1 - Cl 


A similar relation holds for the beam tip 


Sq 2U > = R(£,t)C ni J&ty 


The virtual work of forces and moments at the root and tip nodes must equal the vir- 
tual work caused by generalized forces at the root and tip of the beam 


so that 


f r 

* R 
• ou 

= Q. 

q . 

— 


ii 

H 

m r 

X 1 R 

= Q 

6q 

— 


az 

a2 

H 
Pm 1 

T 

• 6u 

Qa3 < 5 C lci3 

T 

M 

T 

• 6^ 

^a4^a4 


= C TR J Q. 


= R 1 (0,t) J Q 


mM = C TR R T (£,t)^Q 24 


where R T is the transpose of R and the Q's are the generalized forces from the 
discretization of the beam (i.e., coefficients of the 5q f s). 

The element mass matrix may be written explicitly from equation (39) by collect- 
ing coefficients of F u and . After simplification, and ignoring rows and columns 
associated with frame degrees of freedom, the mass matrix can be obtained from 
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^ T 

6u- 


<K V 


60 3 

J 


m6 ij 


3 K 4 

C, . e, . — ff m 

ki kja 3u " a 

p 

3k 3 3k 3 3k^ 3k^ 

Ihl 77 a< + k s^ 7 

a 0 a 3 


symmetric 


m i C 2i “ m 2 C ii 


9k , 


r o 

ii 




•* \~ 


( 70 ) 


To obtain the discretized matrix substitute the shape functions, equation (45), into 
equation (70) and integrate over the element using, for example, Gauss-Legendre quad- 
rature. The dimension of the discretized matrix will depend on the number of internal 
degrees of freedom. Its contribution to the system in terms of nodal degrees of free- 
dom is straightforward and can be left to the reader to determine. Similarly, ele- 
ments of the gyroscopic matrix can be calculated. 


These equations can be programmed either for finite-element computation as they 
are or written in a more explicit matrix form. When the equations are linearized 
about static equilibrium determined from nonlinear static equations, a convenient 
approach is to calculate the mass and gyroscopic matrices explicitly as above and 
solve for the stiffness matrix by numerically perturbing the total static general- 
ized force. 


7 . RESULTS 


In this section, two sets of numerical results are presented along with corre- 
sponding experimental data . The numerical results were obtained by exercising a pre- 
liminary version of a mult ibody /f in ite-element program presently under development in 
which the beam element formulation outlined in section 6 is implemented. The examples 
were set up for calculation using only the properties given in this report. 

The first set of data concerns a cantilevered beam loaded with a tip weight 
(ref. 15). The properties are tabulated in table 1, and the experimental configura- 
tion is described in detail in reference 15. The beam was modeled with one element 
and a sufficiently large number of polynomials to achieve convergence. Results for 
transverse tip deflections u 1 and u 2 and tip rotation 0 3 are shown in figures 4-6 
along with experimental data. The agreement is good, much better than in reference 15. 
This good agreement is achieved, however, without ad hoc modifications of the equa- 
tions based on the values of beam stiffnesses as is done in reference 4. 

The second set of data concerns a cantilevered rotating beam in axial airflow. 

The air flowing down through the rotor plane is induced by the thrust at collective 
pitch setting 0 Q for a two-bladed rotor. Bending and torsion moments were mea- 
sured near the beam root for variable thrust settings. The program was set up to 
calculate air loads based on a quasi-steady aerodynamic formulation by Greenberg 
(ref. 16) and on uniform induced inflow determined from momentum theory. (The 
experimental configuration is described in a forthcoming publication by Sharpe 
(Sharpe , D. L.: An Experimental Investigation of the Flap-Lag-Torsion Aeroelastic 

Stability oi a Small-Scale Hinge less Helicopter Rotor in Hover. NASA TP, to be pub- 
lished in 1985). The rotor properties are tabulated in tables 2 and 3.) The rotor 
blade was modeled with two uniform elastic segments and rigid masses lor the 
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remainder of the structure. Results for a droop angle of 0° and precone angles of 
0° and 5° are shown in figures 7 and 8. The agreement, again, is very good. 


CONCLUSIONS 


Nonlinear beam kinematics for small strains and large rotations have been 
developed and applied to the dynamic analysis of a pretwisted, rotating beam element. 
There are no explicit restrictions on rotation caused by deformation in these equa- 
tions— only the extensional strain of the elastic axis is required to be small rela- 
tive to unity. The only restriction on the magnitudes of the orientation angles 
used in describing the cross-section orientation is that they remain less than 90°. 
For applications of the kinematics where larger rotations may be encountered, a 
method of overcoming the restriction on the magnitude of rotation, which utilizes 
Rodrigues parameters, is presented in the appendix. 

In order to be applicable to all existing rotor/hub configurations in helicop- 
ters, the analysis needs to be incorporated into a hybrid multibody/finite-element 
program. This incorporation is under development. Useful future extensions include 
constitutive equations for composite beams and effects of shear deformation, warping 
rigidity, and initial curvature. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, California 94035, January 18, 1985 
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APPENDIX 


In this appendix the direction cosines and moment strains are expressed in 
terms of Rodrigues parameters. To remove the sign ambiguities of the development in 
the text, the direction cosines should be left in terms of u^, i = 1,2,3 so that 

C 3i * < S 3i + U i )/S ' (A1 > 


This is not necessary in the text because the use of orientation angles limits the 
rotations to be less than 90°; furthermore, it makes the computation slightly more 
involved. For the use of Rodrigues parameters, we make use of similar relationships 
derived in references 8 and 10. The direction cosines in terms of Rodrigues 
parameters c|>^ are 


C. . 


[(1 - ♦*♦*/*) 6 ±j + ^j/2 + e.. k <f» k ]/(l + W 4) 


(A2) 


After much algebraic manipulation, <f> a can be eliminated in terms of the third row 
of C 


<2e a63 C 38 


*3 C so )/d 


+ 


c ) 

3 3 


(A3) 


The use of equation (Al) then yields 

4 > a = (2e ag3 u B + ♦ 3 u a )/( 1 + s ' + U P 


(A4) 


which goes to infinity only when the beam rotations due to deformation reach 
180°. 


The moment strains, as simplified for a straight beam based on those of Reissner 
(ref . 8 ) , are 


K . 
1 


t( hj + e ijkV 2)(f, k ]/(1 + W 4) 


(A5) 


which can be expressed in terms of u|, <J> , u![, and ^3 by differentiation and sub- 
stitution of equation (A4) into equation (A5) . It should be noted that the reference 
basis used by Reissner corresponds to the principal axes of the local undeformed 
beam, instead of the principal axes of the root used herein. Equation (A5) thus 
differs somewhat from Reissner* s equation (49) in reference 8 . 
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TABLE 1.- CANTILEVERED BEAM LOADED WITH TIP WEIGHT 


Property 

Value 

— 

Metric 

Standard 

E 

7.2919X10 10 N/m 2 

(10. 576x10 s lb/in. 2 ) 

G 

3. 022 x 10 1 0 N/m 2 

(4.383x10 s lb/in. 2 ) 

p 

2807 kg/m 3 

(0.1014 lb/in. 3 ) 

ii 


(Ec 3 t/12) 

i 2 


(Ect 3 /12) 

c 

1.270 cm 

(0.4999 in.) 

t 

.3178 cm 

(0.1251 in.) 

L 

i 

50.76 cm 

(19.985 in.) 


TABLE 2.- ROTOR BLADE PROPERTIES (NACA 0012 
AIRFOIL WITH TWO BLADES) 


Rotor diameter, m 1.923 

Blade length (L) , m .870 

Hub offset, %R 9.51 

Chord, cm 8.64 

Taper 0 

Twist 0 


Maximum tip Reynolds number . 600,000 
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TABLE 4.- BLADE MASS PROPERTY DISTRIBUTION 


Inboard 

station, 

r/R 

Outboard 

station, 

r/R 

Mass/length, 

kg/m 

Polar moment of 
inertia/ length , 
kg-m 2 /m 

0.0185 

0.0215 

5.214 

— 

.0215 

.0374 

.0214 

— 

.0374 

.0407 

5.418 

5.827xl0 -3 

.0407 

.0440 

10.010 

7.073x10 3 

.0440 

.0456 

12.745 

4.715xl0 -3 

.0456 

.0555 

9.969 

6.317xl0~ 3 

.0555 

.0608 

5.265 

1.468xl0 -3 

.0608 

.0634 

2.663 

2.082xl0 -3 

.0634 

.0951 

2.429 

2.082x10 3 

.0951 

1.0000 

.343 

2.062xl0 -lt 
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Figure 6.- Geometric twist at tip (deg). 
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